!\Arguments
!  WHICH   Character*2.  (Input)
!          'LM' -> sort X into increasing order of magnitude.
!          'SM' -> sort X into decreasing order of magnitude.
!          'LR' -> sort X with real(X) in increasing algebraic order
!          'SR' -> sort X with real(X) in decreasing algebraic order
!          'LI' -> sort X with imag(X) in increasing algebraic order
!          'SI' -> sort X with imag(X) in decreasing algebraic order
!  N       Integer.  (INPUT)
!          Size of the arrays.
!  X       complex array of length N.  (INPUT/OUTPUT)
!          This is the array to be sorted.
!  Y       complex array of length N.  (INPUT/OUTPUT)
!-----------------------------------------------------------------------
subroutine csortc(which, n, x, y)
    implicit none
    character*2 which
    integer n
    ! Array Arguments
    complex x(0:n-1), y(0:n-1)
    ! Local Scalars
    complex temp
    real temp1, temp2
    integer igap, i, j

    igap = n / 2
    if (which == 'LM') then
        ! Sort X into increasing order of magnitude.
10      continue
        if (igap == 0) return
        do i = igap, n - 1
            j = i - igap
20          continue
            if (j < 0) cycle

            temp1 = abs(x(j))
            temp2 = abs(x(j + igap))
            if (temp1 > temp2) then
                temp = x(j)
                x(j) = x(j + igap)
                x(j + igap) = temp

                temp = y(j)
                y(j) = y(j + igap)
                y(j + igap) = temp
            else
                cycle
            end if
            j = j - igap
            go to 20
        end do
        igap = igap / 2
        go to 10
    else if (which == 'SM') then
        ! Sort X into decreasing order of magnitude.
40      continue
        if (igap == 0) return
        do i = igap, n - 1
            j = i - igap
50          continue
            if (j < 0) cycle

            temp1 = abs(x(j))
            temp2 = abs(x(j + igap))
            if (temp1 < temp2) then
                temp = x(j)
                x(j) = x(j + igap)
                x(j + igap) = temp

                temp = y(j)
                y(j) = y(j + igap)
                y(j + igap) = temp
            else
                cycle
            endif
            j = j - igap
            go to 50
        end do
        igap = igap / 2
        go to 40
    else if (which == 'LR') then
        ! Sort XREAL into increasing order of algebraic.
70      continue
        if (igap == 0) return
        do i = igap, n - 1
            j = i - igap
80          continue
            if (j < 0) cycle

            if (real(x(j)) > real(x(j+igap))) then
                temp = x(j)
                x(j) = x(j + igap)
                x(j + igap) = temp

                temp = y(j)
                y(j) = y(j + igap)
                y(j + igap) = temp
            else
                cycle
            end if
            j = j - igap
            go to 80
        end do
        igap = igap / 2
        go to 70
    else if (which == 'SR') then
        ! Sort XREAL into decreasing order of algebraic.
100     continue
        if (igap == 0) return
        do i = igap, n - 1
            j = i - igap
110         continue
            if (j < 0) cycle
            if (real(x(j)) < real(x(j+igap))) then
                temp = x(j)
                x(j) = x(j + igap)
                x(j + igap) = temp

                temp = y(j)
                y(j) = y(j + igap)
                y(j + igap) = temp
            else
                cycle
            endif
            j = j - igap
            go to 110
        end do
        igap = igap / 2
        go to 100
    else if (which == 'LI') then
        ! Sort XIMAG into increasing algebraic order
130     continue
        if (igap == 0) return
        do i = igap, n - 1
            j = i - igap
140         continue
            if (j < 0) cycle
            if (aimag(x(j)) > aimag(x(j+igap))) then
                temp = x(j)
                x(j) = x(j + igap)
                x(j + igap) = temp

                temp = y(j)
                y(j) = y(j + igap)
                y(j + igap) = temp
            else
                cycle
            endif
            j = j - igap
            go to 140
        end do
        igap = igap / 2
        go to 130
    else if (which == 'SI') then
        ! Sort XIMAG into decreasing algebraic order
160     continue
        if (igap == 0) return
        do i = igap, n - 1
            j = i - igap
170         continue
            if (j < 0) cycle

            if (aimag(x(j)) < aimag(x(j + igap))) then
                temp = x(j)
                x(j) = x(j + igap)
                x(j + igap) = temp

                temp = y(j)
                y(j) = y(j + igap)
                y(j + igap) = temp
            else
                cycle
            end if
            j = j - igap
            go to 170
        end do
        igap = igap / 2
        go to 160
    end if
end
